{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "ee3a16c5-da70-4d74-a087-ef9b5e0801f7",
"metadata": {
"editable": true,
"nbsphinx": "hidden",
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"import numpy as np\n",
"\n",
"from mpi4py import MPI\n",
"\n",
"from triqs.plot.mpl_interface import plt, oplot, oplotr, oploti\n",
"plt.rcParams[\"figure.figsize\"] = (6,4) # set default size for all figures"
]
},
{
"cell_type": "markdown",
"id": "2105a30a-00d6-4a01-894e-31446ea063db",
"metadata": {},
"source": [
"# Semi-infinite chain\n",
"\n",
"This tutorial shows how to setup a minimal impurity model, that of a single fermionic state, and solve it using `triqs_soehyb`. As the first example we pick the analytically solvable case of a semi-infinte chain. (By having the exact solution we will be able to see how the bold hybridization expansion converges as a function of expansion order.)\n",
"\n",
"In the semi-infinite chain with nearest neighbour hopping $t = 1$ the site with only one neighbour has a single particle Green's function with semi-circular density of states (DOS). Thus the chain can be reformulated as an impurity problem of a single site (with zero energy) coupled to a hybridization function with semi-cicular DOS."
]
},
{
"cell_type": "markdown",
"id": "8c39feb5-50e7-4945-8c6f-69b3ccb6cc0b",
"metadata": {},
"source": [
"## Initialization\n",
"\n",
"First we spawn a solver instance and setting up a few general system parameters."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "3bea90ac-4e06-4dc8-965a-2e1d013de004",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Warning: could not identify MPI environment!\n",
" ___ ___ _ ___ _____\n",
" / __| ___| __|__| || \\ \\ / / _ )\n",
" \\__ \\/ _ \\ _|___| __ |\\ V /| _ \\\n",
" |___/\\___/___| |_||_| |_| |___/ [github.com/TRIQS/soehyb]\n",
"\n",
"beta = 5.0, lamb = 1.00E+01, eps = 1.00E-08, N_DLR = 13\n",
"fundamental_operators = [1*c('0',0)]\n",
"H_mat.shape = (2, 2)\n",
"H_loc = 0\n",
"\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting serial run at: 2025-08-17 22:19:09.421767\n"
]
}
],
"source": [
"from triqs_soehyb import Solver \n",
"\n",
"S = Solver(\n",
" beta=5.0, # Inverse temperature\n",
" gf_struct=[('0', 1)], # Green's function structure 1st index: name, 2nd index: dimension of subspace\n",
" eps=1e-8, # Accuracy of Discrete Lehmann Representation (DLR) used for imaginary time response functions\n",
" w_max=2.0, # DLR freqiency cut-off (the spectrum of the model must be in the range [-w_max, +w_max]\n",
" )\n"
]
},
{
"cell_type": "markdown",
"id": "da5e67dc-b514-4010-9cd8-9d580cd74c60",
"metadata": {},
"source": [
"## Hybridization function\n",
"\n",
"In the solver instance the hybridization function is located at `S.Delta_tau` and we set it to have the semi-circular shape of the semi-infinite chain."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "11d8861b-a877-4d03-8092-736728b2f244",
"metadata": {},
"outputs": [],
"source": [
"from triqs.gf import make_gf_dlr_imtime, make_gf_dlr_imfreq, SemiCircular\n",
"\n",
"for bidx, delta_tau in S.Delta_tau:\n",
" delta_w = make_gf_dlr_imfreq(delta_tau)\n",
" delta_w << SemiCircular(2.0)\n",
" delta_tau[:] = make_gf_dlr_imtime(delta_w)"
]
},
{
"cell_type": "markdown",
"id": "cfb335f6-5e15-4bad-b004-e113e17cb073",
"metadata": {},
"source": [
"## Local many-body Hamiltonian\n",
"\n",
"The local Hamiltonian of the single-site is in this case trivial with a single state at zero energy."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "c20ec536-a9ac-4cd5-987c-f1a2e181bee0",
"metadata": {},
"outputs": [],
"source": [
"from triqs.operators import c, c_dag\n",
"\n",
"h_int = 0.0 * c_dag('0',0) * c('0',0)"
]
},
{
"cell_type": "markdown",
"id": "593304cd-ed71-44b9-973b-ba7e15e9a7b9",
"metadata": {},
"source": [
"## Run solver\n",
"\n",
"Finally we run the solver for a few different expansion orders."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "06ca48e3-2414-4e7b-86e5-037b77796c90",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"AdaPol: Hybridization fit accuracy 3.39E-07, using 5 poles.\n",
"(Order, N_Diags) = [(1, 1)]\n",
"max_order = 1\n",
"\n",
" iter | conv | Z-1 \n",
"------+----------+-----------\n",
" 1 | 3.28E-01 | +4.44E-16\n",
" 2 | 8.18E-02 | +1.89E-14\n",
" 3 | 1.96E-02 | +1.91E-14\n",
" 4 | 4.39E-03 | -1.22E-14\n",
" 5 | 9.45E-04 | +3.06E-14\n",
" 6 | 1.98E-04 | -4.44E-15\n",
" 7 | 4.06E-05 | -1.98E-14\n",
" 8 | 8.23E-06 | -4.55E-14\n",
" 9 | 1.65E-06 | -1.28E-14\n",
" 10 | 3.31E-07 | +1.78E-14\n",
" 11 | 6.58E-08 | -1.50E-14\n",
" 12 | 1.31E-08 | -8.88E-16\n",
" 13 | 2.58E-09 | -1.95E-14\n",
" 14 | 5.11E-10 | +3.11E-15\n",
"\n",
"Timing: incl. excl.\n",
"------------------------------------------------------------\n",
"Dyson equation: 0.194 0.194 44.5% |-----------------|\n",
"Hybridization compression (AAA): 0.079 0.079 18.0% |------|\n",
"Pseudo-particle self-energy: 0.006 0.006 1.3% ||\n",
"Single particle Green's function: 0.000 0.000 0.0% |\n",
"Other: 0.158 0.158 36.2% |-------------|\n",
"------------------------------------------------------------\n",
"Total: 0.437 100.0%\n",
"\n",
"AdaPol: Hybridization fit accuracy 3.39E-07, using 5 poles.\n",
"(Order, N_Diags) = [(1, 1), (2, 10)]\n",
"max_order = 2\n",
"\n",
" iter | conv | Z-1 \n",
"------+----------+-----------\n",
" 1 | 3.28E-01 | +4.44E-16\n",
" 2 | 8.18E-02 | +9.10E-15\n",
" 3 | 1.96E-02 | +2.22E-14\n",
" 4 | 4.39E-03 | -6.88E-15\n",
" 5 | 9.45E-04 | +2.89E-15\n",
" 6 | 1.98E-04 | +1.93E-14\n",
" 7 | 4.06E-05 | -9.88E-15\n",
" 8 | 8.23E-06 | -5.55E-15\n",
" 9 | 1.65E-06 | +1.11E-15\n",
" 10 | 3.31E-07 | -1.04E-14\n",
" 11 | 6.58E-08 | +7.55E-15\n",
" 12 | 1.31E-08 | +1.62E-14\n",
" 13 | 2.58E-09 | -4.22E-15\n",
" 14 | 5.11E-10 | -3.67E-14\n",
"\n",
"Timing: incl. excl.\n",
"------------------------------------------------------------\n",
"Dyson equation: 0.380 0.380 42.1% |----------------|\n",
"Hybridization compression (AAA): 0.093 0.093 10.3% |---|\n",
"Pseudo-particle self-energy: 0.170 0.170 18.9% |-------|\n",
"Single particle Green's function: 0.001 0.001 0.1% |\n",
"Other: 0.258 0.258 28.6% |----------|\n",
"------------------------------------------------------------\n",
"Total: 0.902 100.0%\n",
"\n",
"AdaPol: Hybridization fit accuracy 3.39E-07, using 5 poles.\n",
"(Order, N_Diags) = [(1, 1), (2, 10), (3, 100)]\n",
"max_order = 3\n",
"\n",
" iter | conv | Z-1 \n",
"------+----------+-----------\n",
" 1 | 2.95E-01 | +0.00E+00\n",
" 2 | 5.22E-02 | -1.52E-14\n",
" 3 | 6.72E-03 | +2.60E-14\n",
" 4 | 1.07E-03 | +2.62E-14\n",
" 5 | 2.58E-04 | -1.31E-14\n",
" 6 | 6.04E-05 | +9.10E-15\n",
" 7 | 1.14E-05 | -6.88E-15\n",
" 8 | 1.57E-06 | +2.13E-14\n",
" 9 | 2.18E-07 | +1.38E-14\n",
" 10 | 5.39E-08 | +8.88E-16\n",
" 11 | 1.26E-08 | +1.78E-15\n",
" 12 | 2.47E-09 | -2.22E-15\n",
" 13 | 3.63E-10 | +5.55E-15\n",
"\n",
"Timing: incl. excl.\n",
"------------------------------------------------------------\n",
"Dyson equation: 0.793 0.793 18.8% |-------|\n",
"Hybridization compression (AAA): 0.117 0.117 2.8% ||\n",
"Pseudo-particle self-energy: 2.640 2.640 62.6% |------------------------|\n",
"Single particle Green's function: 0.193 0.193 4.6% |-|\n",
"Other: 0.472 0.472 11.2% |---|\n",
"------------------------------------------------------------\n",
"Total: 4.215 100.0%\n",
"\n",
"AdaPol: Hybridization fit accuracy 3.39E-07, using 5 poles.\n",
"(Order, N_Diags) = [(1, 1), (2, 10), (3, 100), (4, 1000)]\n",
"max_order = 4\n",
"\n",
" iter | conv | Z-1 \n",
"------+----------+-----------\n",
" 1 | 3.03E-01 | +0.00E+00\n",
" 2 | 5.92E-02 | -5.77E-15\n",
" 3 | 9.80E-03 | +5.33E-15\n",
" 4 | 1.14E-03 | +8.88E-16\n",
" 5 | 1.69E-04 | -7.44E-15\n",
" 6 | 3.97E-05 | +6.00E-15\n",
" 7 | 9.47E-06 | -2.75E-14\n",
" 8 | 1.90E-06 | -2.31E-14\n",
" 9 | 3.12E-07 | -8.10E-15\n",
" 10 | 3.57E-08 | -1.13E-14\n",
" 11 | 5.58E-09 | +2.00E-14\n",
" 12 | 1.31E-09 | -4.00E-15\n",
" 13 | 3.15E-10 | -1.11E-15\n",
"\n",
"Timing: incl. excl.\n",
"------------------------------------------------------------\n",
"Dyson equation: 1.172 1.172 1.2% |\n",
"Hybridization compression (AAA): 0.162 0.162 0.2% |\n",
"Pseudo-particle self-energy: 88.449 88.449 93.8% |-------------------------------------|\n",
"Single particle Green's function: 3.736 3.736 4.0% |-|\n",
"Other: 0.788 0.788 0.8% |\n",
"------------------------------------------------------------\n",
"Total: 94.308 100.0%\n",
"\n"
]
}
],
"source": [
"from h5 import HDFArchive\n",
"\n",
"max_order = 4\n",
"\n",
"for order in range(1, max_order+1):\n",
" S.solve(h_int=h_int, order=order, maxiter=20)\n",
"\n",
" with HDFArchive(f'data_order_{order}.h5', 'w') as A: \n",
" A['S'] = S"
]
},
{
"cell_type": "markdown",
"id": "aa2f25aa-9993-47ed-8612-58d292a54357",
"metadata": {},
"source": [
"## Visualization\n",
"\n",
"To see the convergence with expansion order we plot the resulting single-particle Green's function."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "ce78e650-b9bb-4b65-8209-55626b9cb583",
"metadata": {},
"outputs": [],
"source": [
"from triqs.plot.mpl_interface import oplot, oplotr, oploti, plt\n",
"\n",
"def plot_dlr_imtime(g_tau, label, n_tau=400, marker='x', linestyle='-', color=None):\n",
"\n",
" from triqs.gf import make_gf_imtime\n",
"\n",
" g_tau_fine = make_gf_imtime(g_tau, n_tau=n_tau)\n",
"\n",
" if color is None:\n",
" color = plt.plot([], [], linestyle+marker, label=label)[0].get_color()\n",
" label = None\n",
"\n",
" oplotr(g_tau, marker=marker, label=None, color=color)\n",
" oplotr(g_tau_fine, label=label, color=color, linestyle=linestyle)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "b2b58b70-083c-4c07-8d7d-5a6da7e261cf",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--> Loading: data_order_1.h5\n",
"--> Loading: data_order_2.h5\n",
"--> Loading: data_order_3.h5\n",
"--> Loading: data_order_4.h5\n"
]
}
],
"source": [
"import glob\n",
"\n",
"filenames = np.sort(glob.glob('data_order_*.h5'))\n",
"\n",
"results = []\n",
"for filename in filenames:\n",
" print(f'--> Loading: {filename}')\n",
" with HDFArchive(filename, 'r') as A:\n",
" results.append(A['S'])"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "e3f1be9b-4b87-461d-aeed-346948beb746",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for r in results:\n",
" plot_dlr_imtime(r.G_tau['0'], label=f'Order {r.order}', linestyle='-' if r.order != 2 else ':')\n",
"\n",
"delta_tau = results[0].Delta_tau['0']\n",
"plot_dlr_imtime(delta_tau, label='Exact', marker='', linestyle='-.', color='black')\n",
"\n",
"plt.ylabel(r'$G(\\tau)$'); plt.xlabel(r'$\\tau$')\n",
"plt.legend(loc='best'); plt.grid(True); plt.ylim(top=0);"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "1f93ea29-7dd2-4c30-9f3f-b12192a6d615",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for r in results:\n",
" diff_tau = delta_tau - r.G_tau['0']\n",
" diff_tau.data[:] = np.abs(diff_tau.data)\n",
" plot_dlr_imtime(diff_tau, label=f'Order {r.order}',\n",
" linestyle='-' if r.order != 2 else ':')\n",
"\n",
"plt.semilogy([], [])\n",
"plt.legend(loc='best')\n",
"\n",
"plt.ylabel(r'max$|G(\\tau) - G_{exact}(\\tau)|$')\n",
"plt.xlabel(r'$\\tau$')\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"id": "cdf81137-b2a1-4d9b-a7c6-3a8cee56b7a3",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"The convergence rate with expansion order is fast and the error at order 4 is about $10^{-3}$."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.17"
}
},
"nbformat": 4,
"nbformat_minor": 5
}